# Load necessary libraries #
library(devtools)
library(SlicerMorphR)
library(rgl)
library(geomorph)
library(readxl)
library(ks)
##############################################
# Load SlicerMorph Data and Preprocessing #
##############################################
# Choose the Log File from Slicer
# (Ensure a line break is added at the end of the .log file before loading)
SM.log.file <- file.choose()
SM.log <- parser(SM.log.file, forceLPS = TRUE)
# Extract Landmark Data
Data <- SM.log$LM
# Create geomorph dataframe
datagdf <- geomorph.data.frame(landmarks = Data)
# Perform Procrustes Superimposition - PCoords stand for Procrustes coordinates
Pcoords <- gpagen(datagdf$landmarks)
# Save Procrustes Coordinates
save(Pcoords, file = "Pcoords.bin")
##############################################
# Semi-Landmark Sliding (SSL vs HL) #
##############################################
# Extract Semi-Landmarks Data
Data2 <- SM.log$semiLMs
# SSL: Sliding using **bending energy minimization**
PcoordsBent <- gpagen(A = Data, surfaces = as.numeric(Data2), ProcD = FALSE, print.progress = TRUE)
save(PcoordsBent, file = "PcoordsBent.bin")
# HL: Sliding using **Procrustes distance minimization**
PcoordsProcDist <- gpagen(A = Data, surfaces = as.numeric(Data2), ProcD = TRUE, print.progress = TRUE)
save(PcoordsProcDist, file = "PcoordsProcDist.bin")
# Extract and Save Centroid Size (Csize)
write.table(PcoordsBending$Csize, file = "PcoordsList_SSL.tsv", sep = "\t")
write.table(PcoordsProcDist$Csize, file = "PcoordsList_HL.tsv", sep = "\t")
##############################################
# Merge with Master Dataset #
##############################################
# Load Pcoords Lists
PcoordsList_SSL <- read.table(file = 'PcoordsList_SSL.tsv', sep = '\t', header = TRUE)
PcoordsList_HL <- read.table(file = 'PcoordsList_HL.tsv', sep = '\t', header = TRUE)
# Load Master Dataset
Masterfile <- read_excel("path/to/YourMasterfile.xlsx")
# Merge Pcoords with Master Dataset
Merge_SSL <- merge(PcoordsList_SSL, Masterfile, by = "LabID", sort = FALSE)
Merge_HL <- merge(PcoordsList_HL, Masterfile, by = "LabID", sort = FALSE)
# Save Merged Tables
write.csv(Table_SSL, file = "TursiopsSpec_SSL.csv", row.names = FALSE)
write.csv(Table_HL, file = "TursiopsSpec_HL.csv", row.names = FALSE)
##############################################
# Create and Save Final Analysis Dataset #
##############################################
# SSL Dataset
spec_SSL <- read.csv("TursiopsSpec_SSL.csv", header = TRUE)
tursiopsBent.dt <- list(gpa.sh = PcoordsBent$coords, cs = PcoordsBent$Csize, spec = spec_SSL)
save(tursiopsBent.dt, file = "tursiopsBent.dt.bin")
# HL Dataset
spec_HL <- read.csv("TursiopsSpec_HL.csv", header = TRUE)
tursiopsProcDist.dt <- list(gpa.sh = PcoordsProcDist$coords, cs = PcoordsProcDist$Csize, spec = spec_HL)
save(tursiopsProcDist.dt, file = "tursiopsProcDist.dt.bin")
##############################################
# Symmetry Analysis (SSL vs HL) #
##############################################
## Define Landmark Pairs for Symmetry Analysis
# SSL: Surface Semi-Landmarks - Replace by your own landmarks
lm.pairs_SSL <- matrix(c(
1:104, 106:109, 111:190, 192:195, 197, 198, 200:217, 219:236,
238:287, 290:299, 301, 302, 304:403, 405:434, 436:445, 448:473,
475, 476, 479:582, 584:621, 624, 625, 627:640, 642:661, 663:760
), ncol = 2, byrow = TRUE)
# HL: Homologous Landmarks - Replace by your own landmarks
lm.pairs_HL <- matrix(c(1:2, 9:28, 31:56, 59:70, 72:73), ncol = 2, byrow = TRUE)
## Perform Symmetry Analysis
# SSL
asym_SSL <- bilat.symmetry(PcoordsBent$coords, ind = dimnames(PcoordsBent$coords)[[3]],
object.sym = TRUE, land.pairs = lm.pairs_SSL, iter = 9)
symm.sh_SSL <- asym_SSL$symm.shape
# HL
asym_HL <- bilat.symmetry(PcoordsProcDist$coords, ind = dimnames(PcoordsProcDist$coords)[[3]],
object.sym = TRUE, land.pairs = lm.pairs_HL, iter = 9)
symm.sh_HL <- asym_HL$symm.shape
# Save Symmetry Data - sym stands for symmetry
YOURNAME_sym_SSL.dt <- list(gpa.sh = PcoordsBent$coords, cs = PcoordsBent$Csize,
symm.sh = symm.sh_SSL, spec = spec_SSL)
save(YOURNAME_sym_SSL.dt, file = "YOURNAME_sym_SSL.dt.bin")
YOURNAME_sym_HL.dt <- list(gpa.sh = PcoordsProcDist$coords, cs = PcoordsProcDist$Csize,
symm.sh = symm.sh_HL, spec = spec_HL)
save(YOURNAME_sym_HL.dt, file = "YOURNAME_sym_HL.dt.bin")
# Load necessary libraries
library(geomorph)
library(readr)
library(ks)
library(rgl)
# Set working directory (modify as needed)
setwd("INSERT_PATH_HERE")
# Load dataset
load("YourName_sym.dt.bin")
######################################################################
# PCA Analysis #
######################################################################
# Perform Principal Component Analysis (PCA)
YourNamePCA <- gm.prcomp(YourName_sym.dt$symm.sh)
# Access eigenvalues
eigenvalues <- YourNamePCA$sdev^2
# Calculate proportion of variance explained by each principal component
variance_explained <- eigenvalues / sum(eigenvalues)
# Print proportion of variance explained by each principal component
print(variance_explained)
######################################################################
# 3D PCA Plot #
######################################################################
# Define colors based on Units
group_colors <- rainbow(length(unique(YourName_sym.dt$spec$YourGroup)))
# 3D scatter plot of first three PCs
plot3d(
YourNamePCA$x[, 1:3],
col = group_colors[as.factor(YourName_sym.dt$spec$YourGroup)],
size = 20
)
# Add legend
legend3d(
"topright",
legend = unique(as.factor(YourName_sym.dt$spec$YourGroup)),
col = group_colors,
pch = 10
)
######################################################################
# 3D Cloud Plot #
######################################################################
# Extract PCA scores
Score <- YourNamePCA$x[, 1:3]
Labels <- factor(x = YourNamePCA$YourGroup)
# Kernel density estimation for group classification
Hscv1.Bd <- Hkda(
x = YourNamePCA$x[, 1:3],
x.group = as.factor(YourName_sym.dt$spec$YourGroup),
bw = "scv",
pre = "sphere"
)
# Perform kernel discriminant analysis
Tt.kda3.Bd <- kda(
x = YourNamePCA$x[, 1:3],
x.group = as.factor(YourName_sym.dt$spec$YourGroup),
Hs = Hscv1.Bd
)
# Plot 3D Cloud
plot(
Tt.kda3.Bd,
colors = c("YourColou1","YourColou2","YourColou3","etc.."),
drawpoints = TRUE,
col.pt = c("YourColou1","YourColou2","YourColou3","etc.."),
box = FALSE,
display = "rgl"
)
######################################################################
# PCA Lollipop Graphs #
######################################################################
PC1 <- plotRefToTarget(
YourNamePCA$shapes$shapes.comp1$min,
YourNamePCA$shapes$shapes.comp1$max,
method = "vector",
mag = 1,
gridPars = gridPar(pt.size = 0.25, pt.bg = "red"))
PC2 <- plotRefToTarget(
YourNamePCA$shapes$shapes.comp2$min,
YourNamePCA$shapes$shapes.comp2$max,
method = "vector",
mag = 1,
gridPars = gridPar(pt.size = 0.25, pt.bg = "blue"))
PC3 <- plotRefToTarget(
YourNamePCA$shapes$shapes.comp3$min,
YourNamePCA$shapes$shapes.comp3$max,
method = "vector",
mag = 1,
gridPars = gridPar(pt.size = 0.25, pt.bg = "green"))
######################################################################
# PERMANOVA Analysis (PAST) #
######################################################################
# Export PCA Scores for further analysis
write.table(YourNamePCA$x, file = "FileName1.tsv", sep = "\t")